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1 Introduction 



Chemotaxis — the response of cells to chemical changes in their surroundings — plays a crucial 
role in the study of organisms. An example of chemotaxis in action is depicted in the remark- 
able video of a neutrophil (white blood cell) chasing a Staphylococcus aureus bacterium (available 



at http://www.biocheinweb.org/neutrophil.shtml). The goal of this project is to develop a model 
and simulation for the movement of the neutrophil in this situation. 

The model we consider is a so-called 'signal centred' approach. Receptors studded throughout the 
neutrophil's cell membrane detect the levels of chemoattractants in the surroundings. This includes 
chemicals secreted by the bacterium nearby. The resulting chemoattractant concentration gradient in 
the neutrophil's surroundings stimulate the cell's movement. Signalling pathways induce protrusions 
in the membrane in the direction of increasing chemoattractant levels; hence towards the bacterium. 
The rest of the cell then follows in that direction through the resulting tension in the cell membrane. 
It should be noted that the neutrophil may respond to the presence of other cells (such as red blood 
cells) as obstructions. However, we shall ignore this factor for the purposes of this model. 

The chemotaxis model of Neilson, Mackenzie, Insall and Webb in [8] acts as our starting point. The 
reaction-diffusion equations for this setup are reproduced below. 

d^b + bVr(t)-v = DbArn)b-nb + rbf adx. (1.1b) 

Jr{t) 

d'c + cVr(t) • V = L>cAr(f)C - rcC + bca. (1-lc) 

where 

• The cell boundary, T{t), is a compact, smooth, connected and oriented curve in for each 
t € [0,T], and moves with velocity v = VfV, where v is the outward normal to T[t). 

• a denotes a local autocatalytic activator (or attractant), with associated decay rate and 
diffusion coefficient Da 

• b denotes a rapidly distributed global inhibitor, with associated decay rate r;, and diffusion 
coefficient 

• c denotes a local inhibitor, with associated decay rate rc and diffusion coefficient Dc 

• In equation (|l.lap . Sa is a saturation coefficient, Sc the Michaelis-Menten constant and ba a basal 
production rate of the activator 

• be in equation (jl.lcp determines the growth of the local inhibitor c in the presence of the activator 
a 

• s incorporates the effect of the external signal and random fluctuations. In [s^ , s is defined as 

s{x,t) = ra[{l + drm^T)) + Rq{1 + drm^V))], (1.2) 

where dr is a positive parameter, RND G (0, 1) is a uniformly distributed random variable, and 
i?o is a non-constant number defined in H related to the local concentration of chemoattractant. 
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In addition, the paper couples this system of reaction-diffusion equations to a PDE controlhng the 
movement of the cell membrane. This PDE is given by 

ut ■ 1^ = Vf — Xk (1-3) 

where u is the velocity field, the outward unit normal to the cell surface T{t), Vf = Kpj-ota with 
-f^prot a positive parameter, k the cortical tension term, and A a term to ensure that the area of the 
cell is controlled: in particular, it is a solution to the non-linear ODE 

dA XoX{A-A, + ^) 

dt- Ao{X + Xo) ^'-^^ 
with Ao and /3 positive constants, A{t) the area of the cell and Aq the initial cell area. 

This project proposes a modified version of the above model by incorporating a new model for the 
neutrophil membrane movement. In Sections [2] and [3l we formulate the mathematical details of the 
normalised chemotaxis model and describe numerical methods to simulate the neutrophil movement. 
Section S] will employ the use of SDEs in establishing escape probabilities of a bacterium being chased 
by a neutrophil. Finally in Section O we consider an empirical model by Li et al. in 0] to compare the 
predictions of our model against experimental observations and consider how the neutrophil's search 
strategy in the absence of a chemoattractant signal improves its chances of finding a bacterium. 



2 Model Formulation 

2.1 Normalisation and Reduction of the Chemotaxis Model 



In Figure 3 of [8[, simulations of the chemotaxis model suggest that, whilst the concentration of local 
activator a varies along the neutrophil cell's membrane, the concentration of global inhibitor h remains 
constant. Furthermore, it appears that the concentration of local inhibitor c follows a similar pattern 
to that of a, but is always a fraction of the value. We want to investigate if the chemotaxis model gives 
rise to this behaviour. Furthermore, these observations suggest that equation (jl.lap alone determines 
the dynamics of the model, and we would like to ascertain whether or not this is the case. Our first 
task is to rescale the reaction-diffusion equations (jl.ip of the Neilson et al. model. By doing this, we 
will be able to reduce the system down to just one key equation. We will then compare our reduced 
system with the data obtained through simulations (as described in 0]). We shall use normalised 
variables x', t', a', h' and c', where 

X = Lx', t = Tt', a = Aa', b = Ab' , and c = Ac' 

and we denote F' be a reparametrisation of F under t' and x'.The normalised form of (jl.ip is then 
given by 

d^,a' + a'Vr'in ■v' = T [j^^r'it')a' + + /j)]; + J^,)^,^) " -'^«') ' (2-la) 

dl,b' + b'Vr'it') ■v' = T (^^Ar,(^t')b' - n {b' - a' dx'^ ^ , (2.1b) 

d'c + cVr'it') ■v' = T ( ^Ar'(t')c' - rcC + bcaj (2.1c) 
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The constants in these equations have specific values in [8| — these are reproduced in Table [TJ We 
would like to choose suitable values for A, L and T to simplify this system so that the dynamics for a', 
b' and c' can be summarised in one, non-linear PDE. Since (|2.1ap appears to be the most complicated, 
involving both a non-linear part and a random term (the value s), we shall try to analyse the other 
two equations in order to see if we can express b' and c' as functions of a' . We can choose: 

• ^ to be the maximum observed value of a, which leaves a' roughly of order 1. According to the 
data of Figure 3 in {sl, we would choose A = 25. 

• L = |r(0)| = Jp^Q^ dx, so that |r'(0)| = 1, and then we can compare cell lengths on a scale basis. 
It is not clear from [8] what L should be since the data provided gives details of the cell shape in 
the later stages of the simulation rather than the initial phase. We guess that a value of L = 0.5 
will do. 



• T to scale with the effect of the Laplace-Beltrami operator in (|2.1ap . Specifically, we would like 
^ = 1, so T = ^ 6.25 X 10^. 

We would like to show the following: 

Theorem 2.1 (Reduced Reaction-Diffusion System). With the above choices of A, L and T , the 
system of equations ()2.ip can be approximately reduced down to 



d:,a' + a'Vp(.) • v' = Ap(.)a' + + +^^1(^.)2,^) " ^Taa' (2.2a) 



6' = f a' dx', (2.2b) 
Jv'{t') 

c' = ^a' ^ 0.385a'. (2.2c) 



We require a few results to prove this. The first of these is the Transport Identity 

^,1 f=[ dlf + fVrit)-v (2.3) 
The second is an approximation result. 

Lemma 2.2. For e > 0, let b{- , t) , be{- , t) : T{t) — t- R be L'^ functions satisfying 

-Ar(t)&, + Cb{be -a) = -e [d^be + beVr(t) ■ v) (2.4a) 
-Ar(t)b + Cb{b-a) = (2.4b) 

Assume that there exist constants Cy,C(i,Ci < oo such that 

• ll^r • ^IIloo < Cv 

• \\dtb + bVr{t)-v\\^^ <Cd 

• |r(t)| < cl, vt > 

Then, for small enough e, be tends to b in L'^ as e — )• 0. 

Proof. We start by subtracting ()2.4bp away from (|2.4ap to get 

-Ar(t) (6, -b) + Chibe -b) = -e [d^b, + b.Vrit) ■ v) ■ (2.5) 
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Now let d := be — b, multiply ()2.5p by d and integrate w.r.t. x over T{t) to obtain 

/ (\Vra)df + Cbd^ + ed{d^b, + beVrit)-v))dx = 0. (2.6) 
Jr{t) ^ ' 

Adding and subtracting ed [d*b + 6Vr(t) • v) gives 

[ (\Vr{t)d\^ + Cbd'^ + <:d{d'd + dVr{t)-v + dtb + bVr(t)-v))dx = 0. (2.7) 
Jr{t) ^ ' 

We then use the Transport Identity ()2.3p with / = ^(i^ to get 



^ I ^ dx + / f |Vrwd|' + + ^Vr(t) • + (9^6 + 6Vr(t) • t;)) dx = 0. (2.8) 



Thus, 



dt ./r(t) 2 dt yr(t) 2 jr(t) 



dx < — I dx + / \^r{t)d 

( ad' + 4 

it) 
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= - /^^^ (^''^^ + -v + ^d [dlb + 6Vr(t) • t^)^ dx 

< {-Cb +^-\\Vvv\\l^) d'dx + e\\dlb + 6Vr(t) -vW^^ j^\d\dx 



- (""^^^I^O / ^^dx + eCsvIrT y Ml^dx 

- ("'^' + I ^"^^ + 2Ci5Cl)) y d^dx 



where we have used Holder's inequality to derive Jp |d| dx < \/\^ Jp \ d\ dx. Now choose e < 



Cb 



. We then obtain 



dt Jr{t) 2 e yr(j) 

By Gronwall's inequality [ill . §1.3], we get 



e\df ^ -Cb f eldP , 
' ' dx < / ^-^dx (2.9) 



ed\ f f ed' 
dx < 



dx ) exp ( ) (2.10) 



r(t) 2 \Jr{o) 2 

Cancel the | factor in ()2.10p and then let e — )• to obtain the result. □ 
We are now ready to show that the reaction-diffusion system can be reduced. 
Proof of theorem \2.1\ Rewrite (|2.1bp into the following form 

dt,b' + b'Vr'tf) ■ v' = J^b^r'anb' -nib' -i a' dx' ) (2.11) 
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where = ^ 5 x 10^ and n = Tn ~ 1-875 x 10^. Upon dividing by Db, we see that the LHS of 
these equations becomes very small in comparison to the rest of the terms, and can thus be taken to 
be zero. We then obtain 

The reasoning behind this is given by Lemma 12.21 — equations ()2.1ip and ()2.12p can be rewritten in 
the form of equations (j2.4ap and (|2.4bp where b' in ()2.1ip is taken to be 6^ in (|2.4ap . C{, = e = i 

and a = ^p^-^^ a dx. 

Now rewrite (j2.12p in variational form to get 

f Vr'(t')b'Vr'(t')V + ^ [ h'v = ^\i a' dx'] [ v (2.13) 
Jr'{t') Db Jr'{t') Db \Jr'{t') J Jr'(t') 

with test function v G H^{r'{t')) =: V. The LHS of (j2.13p is a bounded and coercive bilinear form 
on V and the RHS defines a bounded linear functional on V. We can now apply the Lax-Milgram 
theorem to see that it has a unique solution. Upon choosing 



b' = f a'dx', (2.14) 

we see that it is the unique solution to (|2.13p . Thus it is an approximation of the value of b' under 
the assumptions of Lemma 12.21 



The approach for ()2.1cp is similar. Rewrite the equation as 

d'c' + c'Vr'(i/) • v' - A:Ap/(t/)c' = -f^c' + bcu' (2.15) 

where Dc = -^t" ~ 7, = Trc ~ 8125 and be = Tbc ~ 3125. Upon dividing by f^, we see that 
the LHS of these equations becomes very small in comparison to the rest of the terms, and can thus 
be taken to be zero — although we do not state an explicit result, we hope that a modified form of 
Lemma 12.21 to include a bound on the diffusion term Ap/(^/)c' can justify this approach. We now have 

c' = ha' 0.385a'. (2.16) 

□ 

We have therefore shown that the model does give rise to the behaviour observed in the simulation 
results of and that the system dynamics are contained in the dynamics of a. 

Variational Formulation 

We now derive a weak formulation of (|2.2ap . For notational simplicity, we lose the ' and rewrite (|2.2ap 
into the more compact form 

d*,a + aVr(t) • v = DAr^t^a + /(a) (2.17) 

where /(a) = T ( ^^e)(i+A^(a)^s ) ~ ^"'^ ) ■ Using the Transport Identity (|2.3p . we obtain the follow- 
ing variational formulation: 
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(P^k)- Find a(-,t) E F = H^iGr) such that for almost every t G (0,r), 



4/ a0 + D/ Vr(t)aVr(t)</>= / /(a)(A, (2.18) 

r(t) Jr{t) Jr{t) JT{t) 



dt 



for every G where Gt = U4e[o,T](r(0 x {t}) and /(a) = T ( ^^,^+A,)li+AKay'sa) " 



2.2 Modelling the Neutrophil Membrane Movement 



Recall the neutrophil membrane movement model of [8[ described in Section [TJ One of the problems 
with this model is finding values of A which have to satisfy (|1.4p . a non-linear ODE. Therefore, we 
propose an alternative model that eliminates the Ak term and replaces the formula for Vf with a mean 
curvature flow model, given by 

Vf{x) = -eH{x) + 5a{x) + \. (2.19) 

where e, 5 are small, positive constants, A is a Lagrange multiplier which constrains the area of the 
cell to remain constant and H{x) is the mean curvature at point x, where 

H = Vr(t) • v{x), X G T{t) 

with the outward normal v oriented as shown in Figure [TJ To ensure that the results obtained using 
this model can be compared to those of we choose 6 to be of order Kprot^, in keeping with the 
normalisation of the reaction-diffusion system. The eH term also replaces the cortical tension control, 
represented by the Xk term, in the original mode. We can derive an explicit formula for A. 
Lemma 2.3 (Formula for A). 

A(t) = - 6b (2.20) 

where b is as given in equation ()2.2bp without primes. 
Proof. Since we want the area to remain constant over time, observe that 

= = 1= / ^f= [ -^^(^) + '^"(^) + (2.21) 

a* ctt Jr(t) Jr(t) 

Thus A \T{t)\ = e f-p^^^ H{x)+6 Jp^^^ a{x). However, the Gauss-Bonnet theorem states that Jp^^-^ H{x) = 
2ti when T{t) is a closed curve which is not self-intersecting. Substituting this and ()2.2bp and rear- 
ranging the equation gives the desired result. □ 

Now let X e (M X [0, T],W^) be a parametrisation of T{t). Since = VfU, we get 

X( = —eH{x)v + 5a{x)v + Xv. 

Note that v = (^yx^^ = ]X7' Furthermore, using equation 2.15 of [sl, we can deduce that 



Hu 



-I d / X. 



|X„| dp V|X, 
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We now have 

Xi |Xp| = (^i^ + {5a + A)X^ in [0, 1] x (0, T) (2.22a) 

X(-,0) = Xo in [0,1]. (2.22b) 

We also require that X satisfies the periodicity condition 

X(p,t) = X(p+l,t), p£R, t£ [0,T]. (2.23) 

Suppose that X : M x [0, T] — t- is a smooth solution of ()2.22p and ()2.23p — in particular, 
|Xp| > in [0,1] X [0,r]. Taking the dot product with a test function (p £ ^^^^^([0, 1]; M^) — 
{(j) e {[0,l];M?)\(j){0) = and integrating with respect to p yields the following parametric 

variational form for (|2.22p . 

(P™k)- Given a G H^^^.{[0,1] x [0,r];M), find X G ^^^^^([0,1] x [0,r];M2) such that 

[' [Xt ■ <t>] |Xp| + ^J^I-A dp = r {5a + A) • X^ dp (2.24) 
Jo \-^p\ Jo 

subject to the area of the cell remaining constant, for all (p £ Hp^j,{[0, 1];^?). 

3 Model Numerics 

We now develop a finite element method to numerically solve the reduced system of reaction-diffusion 



equations and simulate the movement of the neutrophil using variational problems (P^^) and (P™^) 
The full non-linear, stochastic system will then be solved to simulate the movement of the neutrophil 
in the absence of chemoattractant (i.e. Rq = 0). 

3.1 Reaction-Diffusion PDE 

3.1.1 Finite Element Approximation 

The smooth, evolving surface T{t) with dT{t) = is approximated by an evolving surface Th{t) with 
drh{t) = 0. Tfi{t) is a polyhedral surface whose vertices {Xj(i)}^^ are taken to sit on r{t) so that 
I'h{t) is an interpolation. Let X'^ : Mx [0, T] — t- be a smooth parametrisation oiTfi{t) with |Xp| > 0, 
and periodicity condition X''{p, t) = X^{p + l,t), < t < T, Vp e R. Then for F : [0, 1] x [0, T] ^ M 
we have 

^-^«^(^'^^-|X,Mp,t)||XpMp,t)r 

This parametrisation will allow us to replace the surface integrals by integrals over the unit interval, 
thus reducing our numerical analysis to a ID finite element method. Let pj = jh, with j = 0, . . . , N , 
be a uniform grid with grid size h = 1/N and define the finite element space 

Vh = {(PeC'^{[0,l];R)\<P\[p^_^^p^^€Pi,j = l,...,N;m=<t>{'^)} (3.1) 
to be the space of piecewise linear, continuous functions. We can now define the semi-discrete problem: 
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(P^). Find a'^{-,t) E Vh such that for ahnost every t G (0,T), 



dt 



X 



dp + D 



1 aj^(/'p 
|X^| 



dp= I f{a 





/in 



dp, 



for every (j){-,t) G Vh- 



Now, denoting the nodal basis functions by {(f)j}jLi, and letting 



TV 



where dim{Vh) = N < oo, we can write the finite element approximation of (j'2.18p as follows. 



dt 



i=i 



dp + DY,Aj 



dp 



X 



dp,i = l,...,N. 



We thus end up with the following system of nonlinear ODEs, 

d 



— ( M(t)a ) + L>S(t)a = M(t)F(a'', i) 



where 



M is the time-dependent mass matrix given by 



S is the time-dependent stiffness matrix 







/ 4>i(t>j 




Jo 





dp, ij = 1, ...,N. 



^0 l^pl 



i,j = 1,...,N. 



a is the time-dependent solution vector 

a(t) = {Ai{t),...,AN{t)), 



F{t) is the time-dependent right-hand side vector 



Jo 



X^ 



dp,i = l,...,N. 



We can solve ()3.4p numerically by discretising it in time via the implicit Euler method, as done 
Let tm = mAt. Then we get 



m+1 -m+l A/r'"Tr"* 



At 



+ S'"+^a'' 



m+l\ m+l 



M'"F'", m = 0,...,M 
M"^ (AtF™ + a"^) 



where A™ := A(tm)- Note that we take F explicitly in time. 



3.1.2 Numerical tests 



Before we attempt to solve (j2.17p numerically, we test this method in a simple setting ; an oscillating 
ellipse of the form 

m := {x = (xi,X2)^ E m2 I +xl = (3.10) 

We consider examples from 

Example 3.1. On the stationary unit circle S"^, the function as{x,t) := e~'^^xiX2 is a solution to 
the surface heat equation dtas — ^s2as = 0. Parametrising the stationary unit circle and deriving a 
variational formulation in terms of this parametrisation as was done in the previous section, we solve 
the resulting problem numerically using the PI finite element method over the domain = [—1,1], 
choosing the final time to be T = 1. Figure [2] shows the numerical solution for a spatial discretisation 
/i = 0.1 and timestep At = 0.0102. Note that for this model problem there is no restriction on our 
choice of timestep as th e sc heme is unconditionally stable. Stability and convergence analysis of this 
scheme can be found in [lC)| . Figure [3] shows the plot of the L? error between the exact solution and its 
finite element approximation at the final timestep t = T for h = 10^^ and At = 10~^. From ^lOi] . the 
implicit Euler finite element approximation of the model problem satisfies the a priori error estimate 

max ||e;rilL2(Q)<C(/i2 + At). (3.11) 

l<m<A/ ^ ' 

By choosing At = 0{h'^), as we have done by our choice of discretisation, we expect quadratic 
convergence in h and this is confirmed by the plot in Figure [3j 

Example 3.2. As another test example, we consider an expanding circle centred at the origin with 
radius r(t) = 0.75 + 5t. The velocity field is given by v{x,t) = ^ and, hence, is purely in the normal 
direction. The function aE{x,t) := exp(g^) ^^^^^p is a solution to dtUE + v ■ Vue — Apag = 0. As 
for the previous test example, the error at the final timestep converges quadratically in h. 



3.1.3 Results 



We now test this numerical method on (|2.1ap by including the source term 

As mentioned earlier in Section I3.1.H we deal with this nonlinear term by taking it explicitly in time. 
The approximations of b and c (respectively ()2.2bp and ()2.2cp ) derived in Section[2]as well as expression 
()1.2p for the stochastic component s(x,t) have been used to compute this source term. The chemical 
constants present in the source term are taken from Table [TJ Given that we have not yet discussed 
numerical methods for the curve evolution, we will restrict ourselves to solving this PDE over the 
stationary circle with radius R = 1. 

Figure H] plots the activator levels at different (rescaled) times, where we have chosen h = and 
At = 400/i^ with final time T = 10. The initial condition for a is chosen to be the 'bump' shaped 
function exp ^— ^ • Now, given that the first term in (|3.12p is of order O (^) whilst the second 

term is of order 0(a), we expect that with our choice of initial condition the activator level would 
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increase rapidly due to the effect of the first term. The second term has the opposite effect when a 
is large enough, by countering the rapid growth induced by the first term. This is exactly the sort of 
behaviour we observe in Figured! 



3.2 Neutrophil Membrane Movement PDE 



3.2.1 Finite Element Approximation 



We shall now apply a finite element method to solve (|2.24|) numerically. Once this is done, this will 
be coupled with the reaction-diffusion PDE discussed previously to simulate the entire system. Define 
the finite element space V/i in the same way as Vh in ()3.ip but with values now in rather than M. 
Motivated by the approach of §4], we derive the corresponding semi-discrete problem satisfied by 
X\ 

(PJf). Given a'^{-,t) £ Vh, find G V/, such that for almost every t £ (0,r), 



H r^n — dn 

X^ 



,5a'' + AM(^- (X^) dp 



for every (p{-,t) G V/j, subject to the area of the cell remaining constant, where A'^ is a discretised 
form of (|2:20D . 

As before let be the standard, piecewise linear, scalar nodal basis. Set 



N 



where Xj(t) G M and let (p = (pje^, k = 1,2, j = 1, . . . , N . Since Y2j=i 4>j = 1) 

N 



Hence, lov i = 1, . . . , N , 

f 

Jo 



5a^ + A'^ ) (/.i ( X: 



Sa'' + A'^ = ^ [6a'' + A^j (t)jek, k = l,2. 



N 1 

•e,dp = ^(<^a'^ + A'^) / 

7 = 1 ■'^ 



Xn 
- > J \ P 



ekdp, k = 1,2. 



We then get 

N .1 



E 







dt 



[Xj] (pj(f)i 



X 



N 



efc dp + J^X, 



i=i 

N 



'^j,p'ri,p 



iV „i 



• ekdp 



efc dp, k = 1,2. 



As before let = mAt, m = 0, . . . , M. The following time discretisation is suggested to derive a 
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numerical scheme for ()2.24p : 




At 



X 



ekdp 



N _ /■! 

,=1 -lo 




h,m 



Note that we take the Lagrangian A'^ exphcitly in time. Now set Xj(f^ 
define, for m = 1, . . . , M, matrices Mx and My by 



efcdp, k = 1,2. 



{xj{tm),yj{tm)), and 



{tr. 



{tr. 



■vh,m 



ei dp, 
e2dp. 



where x"* = (xi(t„t); • • • ^^Nitm))'^ and = (yi(tm), • • • ,yN{'tm))'^ ■ We thus obtain two decoupled 
systems of equations — one for each component of X. 



j^m^m+i ^ £AtS'"y'^+^ = M^y*^ + AtM"^ {6a"' + A"^!) 
where we denote A™' to mean A{tm) and 1 = (1, . . . , l)""". 



(3.13) 
(3.14) 



3.2.2 Numerical tests 



As was done for the reaction-diffusion PDE, we test this numerical method in a simpler setting before 
attempting to solve the original problem. 

Example 3.3. Consider a ball in of radius R{t) in a supercooled liquid of temperature U < 
satisfying the following equation for its surface, 

V = -eH - U 

which is a Gibbs-Thomson law (surface tension —eH with e > 0, and balancing temperature U) 
modified by kinetic undercooling. Then it can be shown that a ball with initial radius Rj shrinks or 
grows depending on whether depending on whether its initial radius is respectively bigger or smaller 
than the critical radius Rc = —jj- Deriving a numerical method for the Gibbs-Thomson law is done 
in a very similar way to our original curve evolution model. 

Figure [5] and Figure [6] show the curve evolution under the Gibbs-Thomson law for e = 1 and U = —1 
(and so Rc = 1) with initial radii Ri = 0.8 and Ri = 1.2 respectively. As expected, the curve shrinks 
for the former and expands for the latter. 



3.2.3 Results 

We now test this numerical method on ()2.19p with a = 0. The full system, including the coupling of 
the curve evolution with the reaction-diffusion PDE, will be dealt with in the next section. Figure [7] 
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plots the curve evolution at different (rescaled) times for an initial curve given by the ellipse 

rO :={x=(xi,X2f GM^ |^+x^ = l} 

where have chosen h = 0.02, At = lOO/i^ and e = 1 with final time T = 10. The initial ellipse should 
theoretically evolve into a circle with the same area. This provides a natural way to test the accuracy 
of the numerical method. We can calculate the area numerically as follows: denote by |^^™| the area 
of the interior of the interpolated curve at timestep m, then we have 

H\ = m= [ 1=/ ^X^'™-z^"''",m = 0,...,M, 

where j/"'™ = "^)^\ corresponding outward normal. The area of the initial curve is given 

explicitly by \^^\ = 27r ~ 6.2832. The numerically computed areas are given as follows: 

\nl\ = 6.2832, \nl\ = 6.2999, | = 6.3117, = 6.2941, 

where the superscript values are the values of m corresponding to the times shown in Figure [71 The 
variation in the computed numerical areas is due to the scheme we chose to consider, which took the 
area constraining Lagrangian term A'^ explicitly in time. 

3.3 P\ill Reaction-Diffusion System 

Having described numerical methods for both the reaction-diffusion PDE and the curve evolution, we 
can now couple the two to simulate the entire system: the autocatalytic attractant ah is computed 
using (j3.9p and then taken explicitly in time in (|3.13p and (j3.14p to evolve Th{t). 

Figure [8] illustrates the progression of the neutrophil under the full system at various times (with 
corresponding activator levels on the right), where we have chosen h = 5 x 10"'^, At = 0{h'^), 
5 = 312.5 and e = 10"^ with final time T = 10. The initial curve is chosen to be the unit circle 
and the initial condition for the activator is given by oq = 20exp (^~o^5SS§^) ('^^ich acts as an initial 
impulse on the left-hand side of the neutrophil) . The model constants are as in Table [TJ 

The result is the formation of a protrusion in the direction of the initial impulse, driving the neutrophil 
to move in the same direction. However, as the figure shows, the advection of the neutrophil results 
in a progressive clustering of the points which will lead to an ill-conditioned problem. Our solver will 
thus have to include a remeshing step. 

3.3.1 Remeshing 

The procedure for mapping Xj to the remeshed point is given as follows: 

X- = + -^^+1 _^ (("Xi - Xi_i) • nj)ni, 1 = 1,... ,N, 

where rij is the outward normal of the segment [Xj_i,Xj]. Figure [9] illustrates this procedure. It is 
important to note that this remeshing step is area-preserving, which we require to ensure that the cell 
area remains constant. 
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This process alone creates an artificial, tangential transport component v!^ that has to be counter- 
balanced in the activator PDE with an advection term in the opposite direction, given by the matrix 

/■I x' 
= / At{v')\p) ■ j-^ Xi,p Xj dp, i,j = 1,...,N 
Jo \-^p\ 

with [v')^ = + ) where is the normal velocity of T}^. We are thus required to solve the new 
system 

This allows us to adapt the activator values accordingly when the remeshing is performed. Figure [TOl 
shows the effects of the remeshing on the progression of the neutrophil under the full system at 
various times (with corresponding activator levels on the right), where we have chosen the same 
discretisation/parameter values and initial conditions as before. 

3.3.2 Remarks 

The results based on the original model of suggest that a "parent" pseudopod would split to give 
rise to two new "child" pseudopods, a process that can be observed in the activator profile. Our results 
do not appear to show any such splitting: the one pseudopod that was created by the initial impulse 
does not result in any further pseudopods being formed. 

The reason for this may in fact be due to our reduction of the original model; in particular, neglecting 
the diffusion term in (jl.lcp . The Diffusion coefficient Dc is greater than Da, indicating that the local 
inhibitor may have an effect over a slightly larger part of the pseudopod than the local activator. We 
argue that this is because the activator's influence is over the tip of the pseudopod, whereas the local 
inhibitor affects both the pseudopod tip and the parts of the membrane connecting the pseudopod 
to the rest of the neutrophil. This suggests that if we want to observe pseudopod splitting then we 
should adapt our model to use ()2.1cp rather than (|2.2cp . 



3.3.3 Extension 

Motivated by the remarks above, the original local inhibitor equation (|2.1cp is incorporated back into 
our model. Figure [TT] shows the resulting simulation of the neutrophil at various times (with corre- 
sponding activator and local inhibitor levels on the right), where we have chosen h = A x 10~^, At = 
0{h^), 5 = 312.5 and e = W'^ with final time T = 10. The initial curve is chosen to be the unit 

0002 ) • '^^^ pseudopod 

splitting discussed in the previous section is now clearly visible, confirming our intuition that the local 
inhibitor c is crucial for this phenomena to arise. 
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4 Probabalistic Approach to Bacterium and Neutrophil Movement 



We now address the effects of cell movement and chemotaxis on the probability of a neutrophil catching 
and neutralising a pathogenic bacterium within the bloodstream. We initially model this game of 'cat- 
and-mouse' on an infinite plane before considering the movement of the neutrophil with respect to 
the chemotaxical effects of the bacterium and then restricting the problem to within a capillary with 
negligible net blood flow. Finally, we consider the same capillary problem, but with an additional 
obstacle (perhaps a neighbouring blood cell) which either engulfs or deflects the bacterium. 

It seems reasonable to model the bacterium's movement as a Brownian Motion since we assume 
it is suspended in a fluid (blood plasma). Thus, pendent in this dynamic liquid, it is consistently 
undergoing collisions with the molecules of that liquid and if the number of collisions is large enough 
then the Central Limit Theorem gives rise to a Brownian Motion. We give this Brownian Motion a 
variance parameter, a, which will depend on the dynamics of the fluid and the mass and dynamics of 
the bacterium itself. 

However, when modelling the movement of the neutrophil, the aforementioned chemotaxical effects 
allow us to be far more deterministic in our approach. For our purposes we first assume that the 
neutrophil is centred at the origin, and any of its movement is imparted onto the movement of the 
bacterium — i.e. we observe the process from the neutrophil's perspective. For example, in two- 
dimensions, if we suppose that the neutrophil undergoes movement with constant velocity and moves 
with speed and direction (1, l)""^ (for every unit of time) then, from its perspective, this is equivalent to 
an inclusion of a drift term in the movement of the bacterium equal to (—1, — l)"*"- Thus, in combination 
with the motion of the lone bacterium given above, the stochastic differential equation (SDE) that 
represents the path of the bacterium in our new frame is now given by 



Owing to the fact that chemotaxis is a process based on the detection of a concentration of particles, 
we suggest that it is a diffusive process; thus its effects decay exponentially as one moves away from 
the source. Therefore, in two dimensions, we suggest the following model equation for the movement 
of the neutrophil with respect to the position of the bacterium, {xi,X2)'^: 



where 9 is some scale factor dependent on the cell's detection of chemoattractants and its subsequent 
movement in that direction. This suggests that the SDE describing the path of the bacterium from 
the perspective of the neutrophil is given by 




or, more generally. 



dXt = hdt + crdBt. 
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In Sections 12.21 and 13.21 we considered how the neutrophil membrane varies over time in response to 
the bacterium's presence with a view to couphng the equations derived here with the reaction-diffusion 
and membrane movement PDEs. For this section, however, we work with the unhkely, but hopefuhy 
not too inaccurate, assumption that the neutrophil membrane maintains a perfectly circular shape 
and is unaffected by the movement of the neutrophil or the dynamics of the surrounding fluid. 

We now turn our attention to the various domains that may be considered. Initially, we consider a 
bacterium trying to escape a neutrophil on an infinite, two-dimensional plane. Using R, the radius of 
the neutrophil, as a convenient scale factor, we assume that the bacterium succeeds in its escape if it 
reaches a distance of kR from the neutrophil, for some k > 1. This provides us with a domain in the 
shape of an annulus, 

n = {x e M.'^IR <x< kR}. 

From here, we can then extend this to a bacterium within a capillary. Finally, we shall attempt to 
add obstacles such as other blood cells within the capillary. 

As mentioned previously, the path of the bacterium can be given by an SDE with a drift term b{Xt,t) 
(examples of which were given earlier) and a constant variance, a. Using the annulus as our domain 
with kR as the escape radius, we propose a method for computing the probability that a bacterium 
is engulfed by the neutrophil given that the bacterium starts at a point /3 G $7°. 

Suppose we also have a function f{x,t) £ ([0,oo) x M), which describes the probability that__the 
bacterium will escape given that it is in position x at time t. Consequently, using Ito's Lemma 



dfiXt, t) = dtfiXt, t)dt + vfiXt, t) ■ { biXt,t)dt + 

= Af{Xt,t)dt + VfiXt,t)- I " \dWt 





dWt 



+ ^^'A/dt 



(4.1) 



where Af(Xt,t) is the generator associated with our SDE [^]. Integrating ()4.ip from to T and taking 
expectations gives 



E[f{Xt,t)]-f{Xo,0)=E / Af{Xs,s)ds 

Jo 

By the martingale property of the Ito Integral, 

E[f{Xt,t)] = f{Xo,0)+E 



+ E 


[I 







^ Vf{Xt,t)- llj<iWt 



(4.2) 



AfiXs,s)ds 



Equivalently, by the Optional Stopping Theorem, for stopping time r = inf{t > 0\Xt ^ fi} (the time 
that the bacterium exits our domain), 

E[/(X„r)] = /(/3,0)+E r Af{Xs,s)ds 

Jo 

where f3 is the starting point of the bacterium. However, we also note that 

E [f{Xr,T)] = F{Xr = R)f{R, t) + F{Xr = kR)f{kR, r). 
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Thus, 

¥{Xr = R)f{R, r) + ¥{Xr = kR)f{kR, r) = 0) + E Af{Xs, s)ds 

Jo 

for every f{x, t) G ([0, oo) x M). So if we can solve the partial differential equation (PDE) Af = 
with boundary conditions f{X) = 1, when |X| = kR, and f{X) = 0, when \X\ = R, this provides us 
with 

V{Xr = kR) = f{P,0) (4.3) 
since f{kR, r) = E [£ Af{Xs, s)ds] = 0. 

However, when solving Af = for a parabolic generator, we need a third, time-related boundary 
condition. To this point, for a time dependent drift term b{Xt,t) we have the following problem: 

dtf + b{x,t)-Vf + ^a''Af = 

/ = in dBniO) x (0, oo) (4-4) 
f = l in dBkRiO) X {0,oo) 

For the third boundary condition, we make the assumption that the time dependence of the drift term 
diminishes over time — i.e. the neutrophil's strategy for chasing the bacterium still depends on the 
position of the bacterium, but becomes less volatile over time. Equivalently, 

b{Xt,t) d{Xt) as oo. 

for some function d that depends only on Xt. Furthermore, we suppose that there exists a time T < oo 
where 

Vi>r, b{Xt,t) = d{Xt), 

Thus our third boundary condition for our parabolic problem is the solution, u{x), to the elliptic 
problem 

d{x) ■ Vu + ^c^^An = 

u = in dBniO) X {0,oo) (4-^) 
u = 1 in dBkRiO) x (0, oo) 
Finally, letting s = —t to obtain a forward heat equation, we derive the parabolic PDE 

dJ-bix,s)-Vf-^a^Af = 

f = in dBniO) X {-T,0) 
f = l in dBkdO) X {-T,0) 
f = u{x) in nx {-T} 

The completion of this formulation of the problem involves showing that the values of solution to this 
PDE, f{x,t), take values between and 1, as required of a probability function, since our solution / 
represents the probability of the bacterium escaping. 
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Theorem 4.1. If f{x,t) satisfies ()4.6p with the given boundary conditions on the domain Q, which 
consists of an annulus with inner radius R and outer radius kR, then Vx G G [— T, 0], 

0</(x,t)<l 

Proof. The partial differential equation in question is equivalent to 

dtf + Lf = Q where L/ = t) • V/ - ^a^A/ 
Now, since dtf + Lf < 0, by the Weak Maximum Principle Q], 

max / = max / = 1 and min / = min / = 0. 

ax[-T,0] 9Cx[-T,0] nx[-T,0] Snx[-T,0] 

Thus, Vx G G [-r,0], < < 1. □ 
With this general framework in place we attempt to find a solution for several special cases. 

4.1 No Neutrophil Movement 

The first and simplest case that we consider is the case on the two-dimensional annulus when there 
are no chemotaxical effects and therefore negligible neutrophil movement. 

Theorem 4.2. For a stationary neutrophil and a bacterium that is modelled by a Brownian Motion 
with variance 1 beginning at a point (5 such that R < |/3| < kR, where R and kR are the cell radius and 
escape radius respectively, and that is almost surely engulfed if it makes contact with the neutrophil, 

1 //3 

F [Bacterium engulfed) = 1 — - — - In ( — 

111 rv \ Xt 

Proof. In this case the bacterium follows the path 

dZt = dBt. 

So the drift term is zero and there is no time dependence. By (j4.3p . we observe that the solution 
to this problem is given by the solution f{x,t) to ()4.4p with no time dependence and the drift term 
b{x, t) equal to zero. Therefore, the probability of escape is given by the solution to the PDE: 

-A/ = 

f = in dBniO) 
f = l in OBkRiO) 

Since we are on the annulus it will be easier to change to radial coordinates and we note that we 
choose a radially symmetric solution. 

r dr \ dr 



Solving this ODE gives 



This gives the probability that the bacterium will escape given that it starts at r. Hence the probability 
that the bacteria is engulfed is given by 1 — /(/?). □ 
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Remark. We also note that since Brownian Motion is recurrent (in the sense of 0, §4]), to a ball 
with area greater than zero in two dimensions, as we increase the escape radius kR we also increase 
the probability that the bacterium is caught. In fact, as A; — )• oo, P(Bacterium engulfed) — )■ 1. With 
infinite time and no obstructions the neutrophil will always engulf the bacterium. 

Similarly, in three dimensions: 

Theorem 4.3. For a stationary neutrophil and a bacterium that is modelled by a Brownian Motion 
in three dimensions beginning at a point (3 such that R < \/3\ < kR where R and kR are the cell 
radius and escape radius respectively, and that is almost surely engulfed if it makes contact with the 
neutrophil, 

f R ' 

^'{Bacterium engulfed) = 1 — I 1 



Proof. In the same way as the previous theorem 



I d ( ^df 



^ r^ dr \ dr 



Solving this ODE gives 




Again, the probability of being engulfed is 1 — /(/?). □ 



4.2 Movement of the Neutrophil 



Now that we have obtained a rough outline of the probability function, we endeavour to find solutions 
in the more general case when 6 7^ (movement of the neutrophil) . Analytic solutions to the equations 
(j4.4p are not obvious and none of the immediate methods for solving PDEs, such as the separation 
of variables technique, produce a solution. Therefore, we move to calculating these solutions com- 
putationally using linear finite element methods and linearly approximated domains on the software 
package DUNE. In Figures [T^ to [T71 we model the probability of a bacterium escaping a neutrophil, in 
a variety of situations, as a function of the starting position of the bacterium relative to the neutrophil 
(centred at the origin of each figure), i.e. a value close to one at a point in these figures implies that 
a bacterium, which starts at that point, is very likely to escape the neutrophil. Throughout these 
figures, the left image displays the PDE solution of the probability function that models the probability 
of a bacterium escaping and the right image displays the 'iso-probability lines' joining possible starting 
points of the bacterium that give an equal probability of the bacterium escaping. 

Figures [T2l 1131 andll4l demonstrate how the choice of drift term affects the bacterium's chances of escape 
on an annulus domain and the PDE representing this, together with its associated computation, is 
discussed in detail in the captions below these figures. The unlikeliness of this particular environment 
is ameliorated by considering the more realistic setting of a capillary in Figure [T5l This case is 
considered by solving equation (|4.4p on a rectangular domain with Dirichlet boundary conditions, 
which are described, together with the domain approximations, more thoroughly in the figure caption. 
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However, in figure [T5l we have assumed that a bacterium can pass through capihary walls. We attempt 
to rectify this by using reflective boundaries to represent an impassable capillary wall in Figure [TU] (this 
too is discussed further in the associated caption) . Finally, in Figure [171 we introduce an obstacle, 
which we suppose represents a passing blood cell. For this scenario, we assume that the bacterium 
can either enter or feed off the cell causing irreparable damage in both cases. Thus, we assume that 
if the bacterium reaches the cell the neutrophil has failed in its attempts to neutralise it. Again, for 
further discussion on how this model was formulated and computed the relevant caption should be 
consulted. 

5 Empirical Model for Neutrophil Motion 

5.1 Experimental Data 

Empirical data from [7] suggest that in the absence of the chemoattractant signal from a bacterium, 
neutrophils move in a roughly zig-zag manner. It is suggested that this could be a strategy for 
improving the efficiency with which the cell searches an area of space for bacteria. We now compare 
the PDE models developed earlier to the experimental data using a model based on its statistical 
analysis, and postulate some ways of determining how efficient the search strategy is compared to, 
say, a simple random walk. The data suggest that the cell moves in an essentially straight line for 
a random time with an Exponential(0.67) distribution, measured in seconds. After each such time 
period, the cell makes a turn with magnitude determined by another Exponential(0.67) distribution 
(in radians). The history of left/right turns can be modelled as a discrete Markov chain in which 
the probability of turning the opposite direction to previously (that is, a left turn followed by a right 
turn or vice versa) is approximately 2.1 times the probability of turning in the same direction as 
before: the analysis of 4822 turns from 12 observed trajectories reveals that opposite-type turning 
pairs outnumbered same-type pairs 3623 to 1559. The cell is assumed to move at a constant speed of 
7.46 X 10-^ms-i Q- 

5.2 Comparison with PDE Model 

We should expect our reaction-diffusion model for chemical activators on the membrane to produce 
results similar to those of this second, less complicated model, backed up by experimental data. It 
is suggested by Li et al. that the cell's movement is controlled by the formation of pseudopods, 
which cause the cell to move in the direction they are extended. The zig-zag path indicates that 
the production of pseudopods should therefore occur in an alternating left-right sequence and this is 
behaviour we can examine in the results of the PDE model. 

The PDE model shows that the formation of a pseudopod A is followed by a period in which the 
pseudopod splits in two. The stochastic forcing term in the reaction-diffusion equation (jl.lap causes 
one of these peaks to eventually dominate and grow into a new pseudopod. This is most likely to 
be the peak closest to ^'s parent. New pseudopods therefore form more often between the two most 
recent extensions, and this accounts for the observed alternating turn behaviour. We postulate that 
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this final aspect of the pseudopods' behaviour may be predicted by the PDE model: see, for example, 
the extensions section. 

5.3 Analysis of the Search Strategy 

The cell's aim should be to quickly locate the chemoattractant signal of a bacterium, so we expect it 
to have evolved a search strategy for doing this in an efficient way. A reasonable requirement is for 
the cell to explore the space rapidly without covering the same region too often or getting stuck in 
one place. For comparison we may consider how the strategy improves upon a typical random walk 
without directional bias. 

An example cell path is given by Figure [TE[ the turn history for the same path is illustrated in 
Figure [T9j The type of motion modelled here is a persistent random walk, since the direction of 
movement is dominated by a slowly-changing parameter 9 around which the frequent turns impose a 
highly variable (p term. We can observe this tendency of the cell to keep moving in roughly the same 
direction by conditioning its motion to begin in the positive x direction and observing the points at 
which it exits a circle of a given radius r (centred at the origin). Figure [20] plots the exit points for 200 
choices of radius at regular intervals from 0.1 to 20. The simulation for each radius was run 100 times 
and all exit points correspond to independent paths. The high concentration of points surrounding 
the positive x-axis clearly indicates the tendency of the cell to keep moving in the direction it started 
on average; however we have also seen how the high frequency of turns within this overall trajectory 
causes the zig-zag motion mentioned earlier. The result is a 'sweeping' procedure across the plane 
that may help the cell to explore a large area efficiently and without going back on itself too often. 

Histograms for the cases r = 20, 50, 100, 200, as shown in Figure [21] (100,000 simulations each) further 
press the point, and we note for these higher radii that the influence of the starting direction on the 
exit point diminishes as the radius increases. This shows that while the cell trajectory has a high level 
of local persistence, globally it does not conflne itself to a particular area of the plane. We conclude 
that it is unlikely for the cell to continue moving in the same average direction for an excessive amount 
of time; this is an important consideration in the real-world application of the model where the space 
is not inflnite and therefore not all areas are equally worth exploring. We see that even though the 
cell in our model is not conflned to a finite space, its tendency is not to continue moving in the same 
direction for an excessive time period. 

Figure [52] shows the mean exit times for each integer radius 10-200. The curve is well-fitted by a 
quadratic y = 0.00328x^ + 0.117x + 0.748. Long-term behaviour is therefore that movement from the 
origin occurs on the timescale r^, like a random walk, whereas short-term behaviour is on the much 
faster timescale r. This further backs the conclusion that the cell explores local regions quickly, but 
dwells for a longer time in larger areas. 

6 Conclusion &; Further Work 

In this report, we have taken the model of [sl and reduced it to a single reaction-diffusion PDE based 
on simulation data provided in the paper. In addition we have also proposed a new model for the 
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movement of the neutrophil membrane using a mean curvature approach. Numerical methods for 
the reaction-diffusion PDE and the neutrophil movement have been developed and implemented to 
simulate the system. These simulations do not support the observations made in 0|, namely the 
pseudopod splitting, suggesting that our reduced model was oversimplified. On reintroduction of the 
local inhibitor PDE, resulting simulations show this pseudopod splitting behaviour. Therefore, our 
model gives the same results as in Q], with the additional benefit of having a simpler model for the 
neutrophil movement, which does not involve solving a non-linear ODE. 

In addition to the PDE approach, we also modelled the neutrophil in its pursuit of bacteria through 
the use of SDEs. On an infinite plane it was suggested that movement along a constant vector would 
be the most efficient strategy for the neutrophil, but on the restricted environments commonly found in 
practice our model has found this particular strategy to be ineffective. Moving on from this conclusion, 
we then found models that included chemotaxical effects to be far more fruitful for the neutrophil in 
these situations, as one would suspect. However, the model showed a neutrophil solely based on 
a response to chemoattractants to be lacking over longer distances; suggesting that a compromise 
between these two strategies would be an optimum evolutionary choice for a neutrophil that may find 
itself in a variety of surroundings and would need to perform accordingly. Results, at least partially, 
supported by our other approaches in the relevant sections. Finally, our model portrayed a very 
different picture when obstacles were added. Though our model is a very rough approximation, this 
perhaps suggests that the neutrophil will be at a slight disadvantage in densely populated areas of 
cells. 

Finally, our simulations based on the experimental data indicated that the cell's search strategy is to 
zig-zag frequently back and forth over an area while travelling in one dominant direction for longer 
periods of time. The high persistence of the cell's motion means that it explores its local area much 
more quickly than a standard random walk while dwelling for longer in larger areas, as we saw from 
the analysis of exit times and exit points for circles of varying radii. The PDE model reflects some of 
the observations regarding pseudopod behaviour. 

One way to extend the work done in this project is to consider an alternative model for the neutrophil 
membrane movement. The natural extension would be to consider a Willmore flow model since this 
minimises an energy functional related to the curvature of the membrane. We can also incorporate 
the bacterium via the full expression of the stochastic term {Rq ^ 0), and the SDE model developed 
in Section S] can be used as a benchmark for comparison. In particular, it would be interesting to 
consider the effect (caused by morphing of the membrane) of chemoattractants (6) on the neutrophil 
movement. Finally, we could compare simulated paths produced by both the PDE model and the 
empirical model using statistical comparison methods. Properties such as the expected exit times 
and exit position probabilities could be used for these comparisons. The PDE model could also be 
extended to better reflect the observed positioning of new pseudopods. 
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Figure 1: Orientation of curve parametrisation and normal vector. 
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Figure 3: error at t = T for Model Problem. 
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Activator levels at various times 
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Figure 4: Local activator levels in the simulation of the reaction-diffusion equations 
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Figure 7: Curve evolution under ()2.19p given a = 0. 
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Activator levels at different times. 



Curve evolution without remeshing. 
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Figure 8: Membrane progression and local activator concentration levels simulated by the full system. 
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Figure 9: Illustration of the remeshing procedure. 
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Activator levels at different times. 



Curve, evolution with remeshing. 
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Figure 10: Membrane progression and local activator concentration levels simulated by the full system 
with remeshing. 
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Activator and local inhibitor levels at t=5*10"- 
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Activator and local Inhibitor levels at t=1.2*10^. 
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Figure 11: Membrane progression and local activator and inhibitor concentration levels simulated by 
the full system including the original equation for local inhibitor c. 
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(a) 



(b) 



Figure 12: We begin with the case when b = (5, 5)"^ and a = 1. This is an example of a neutrophil 
moving fairly quickly in the (1,1)"'", direction unaffected by chemotaxis. In an infinite domain it 
could be argued that this strategy (moving in a constant direction) would be the optimum strategy 
that could be employed by the neutrophil since this strategy would cover the most area in a limited 
amount of time. However, from the figures, on this restricted domain, it appears obvious that a 
bacterium is very likely to escape unless it starts in the direct path of the neutrophil or extremely 
close it. Despite the strong initial argument for this strategy it seems to be largely ineffective on 
more restricted domains. The figure itself shows a solution to a steady state case of a PDE derived 
previously (equation (j4.6p ) to provide the probability of a bacterium escaping on an annulus. Thus, 
it is given by: (5,5)"^ • f{x) + ^A/(x) = in 0; / = in (9^0.2(0) and / = 1 in dBi{0), where 
n = {x£ R'^\R <x< kR}. 
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(a) 



(b) 



Figure 13: For further comparison, We also observe the case when b = (10, 10)"^ and a = I. This 
is an example of a neutrophil moving very quickly in the (1, l)""" direction unaffected by chemotaxis. 
By the original argument for a neutrophil's strategy on an infinite domain this strategy should prove 
even more beneficial. However, on this restricted domain it only serves to strengthen the arguments 
against this strategy. The figure shows the solution to the same PDF as before, derived in Sectiord] 
(equation ()4.6p '). but with the noted change to the diffusion term: (10, 10)""" • f{x) + |A/(x) = in fi; 
/ = in 9^0.2(0) and / = 1 in dBi{0), where Q = {x e M^ji? < x < kR}. 
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(a) 



(b) 



Figure 14: We now consider the chemotaxical effects with the aforementioned suggestion of a drift 
term that induces these effects. We use the drift function suggested in Section [J] with 9 = 0.5: 
b{x) = (0.5a;ie~^l^l~'^\ 0.5a;2e~^l^l~'^^)"^. It should be noted that this function has no empirical or 
theoretical background and was only chosen because we expect this function to confer several of the 
properties we expect from chemotaxis: importantly, a greater effect as the bacterium gets closer to the 
cell; this effect should decay exponentially because it is based on the diffusion of particles. The most 
prominent thing we note is that chemotaxis is a process based upon concentrations of amino acids 
in a fluid and the diffusions of these concentrations so the effects of this process decay exponentially 
with distance. Due to this, in our model at least, the bacterium significantly increases its chances of 
escaping by starting only a relatively small distance further from the neutrophil. In addition to this 
we also note that despite the very limited movement we have put on the cell, with the 9 = 0.5, term 
this strategy seems at least as good as the previous strategies in this environment and requires far less 
movement and therefore energy on the part of the neutrophil. The PDE governing the probability 
distribution here is given by: {0.5xie~^^^^~'^\0.5x2e~^^^^^^^)'^ ■ f{x) + ^A/(x) = in fi, / = in 
a5o.2(0) and / = 1 in a5i(0), where 9. = {x e M^ji? < x < kR]. 
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(a) (b) 

Figure 15: We again consider the chemotaxical effects with the same drift term as Figure O How- 
ever, here we consider the domain of a capihary, which we assume, for the moment, to be per- 
fectly rectangular. Furthermore, we suppose that if the bacterium hits either capillary wall, or 
as before, it reaches an arbitrary distance away from the neutrophil, then it has escaped. This 
escaping can be viewed as a bacterium passing through the capillary wall. The results are simi- 
lar to before, but in this case the relationship between probability of escape and distance is more 
pronounced on this longer domain. As we would expect, the bacterium also stands a greater 
chance of escape if it moves horizontally away from the neutrophil than vertically through the 
same distance. The equations for this situation are as given for Figure [m but instead, / = 1 
in r = {x e R'^\\x\ = VW,-l < xi < 1} U {x £ R'^lxi = -1} U {x e R'^\xi = 1} and 
n = {x e R'^\R <x< kR}. 
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(a) (b) 

Figure 16: We now consider the same situation as Figure [T5l but in this case we assume that the 
capihary wall is a barrier that the bacterium cannot pass through. We model this as a reflective 
boundary. This is done through the use of Neumann boundary conditions on the reflective boundaries 
as in [l[ in addition to the already present Dirichlet conditions. In real-life we would expect the 
boundary to have attributes of both this figure and the previous figure. However, of the two we 
assume this to be the more likely and thus we continue with this approximation whilst keeping in mind 
the results of the previous figure as a distinct possibility that cannot be excluded in any extensive 
future modelling. In this figure, we also add the assumption that the bacterium will take up a larger 
proportion of the capillary and therefore, we limit its horizontal movement even further, but increase 
its size to compensate. The results are similar to before, but in this case horizontal movement of the 
bacterium away from the neutrophil is rewarded less than vertical movement in contradiction to the 
previous figure. The equations governing the escape probabilities in this case are: (O, 0.5x2e~^l^l~^^) • 
fix) + ^Af{x) = in f^; / = in 9^0.5(0); / = 1 in Fi and V • / = in F2, where Ti = {x e R'^\ \x\ = 
VTO, -1 < 2;i < 1}; F2 = {x e M^i = -1} U {x e M?\xi = 1} and Q = {x e R'^\R <x< kR}. 
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(b) 



Figure 17: This time we remove the drift term and, as before, we attempt to model neutrophil and 
bacterium movements in the domain of a capillary with impassable boundary walls; modelled as re- 
flective boundaries. However, on this occasion we introduce an obstacle, which we view as equivalent 
to a passing red blood cell. We assume that the bacterium can enter or feed off the cell causing 
irreparable damage and thus we assume that if the bacterium reaches the cell then the neutrophil 
has failed in its attempts to neutralise it. It is perhaps surprising that the obstacle has such a great 
effect on the neutrophil's performance. In fact, it severely limits the neutrophil's ability to capture 
bacteria in the vicinity of the obstacle and, with several obstacles in place, one would expect the 
neutrophil's effectiveness to be drastically reduced. In this case we change both the domain and the 
drift function from the equation used in previous figures, as well as those derived in Section HI We 
assume that the capillary is perfectly rectangular, as in previous figures, and we represent the ob- 
stacle as a square removed from the domain. Thus, the equations characterising this situation are: 
iA/(x) = in ri; / = in 3So.5(0); / = 1 in Ti; / = 1 in Tobstade and V • / = in r2 where 

Ti = {X G M2| = ^,-1 < XI < 1}; T2 = {X £ R^\xi = -1} U {x £ M^l = 1}; ^obstacle = 

the perimeter of the square spanned by the points (0.8, 0.8); (0.9, 0.8); (0.8, 0.7) and (0.9, 0.7) and 
n = {xe R'^\R <x< kR} 
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Simulation of a typical neutrophil path in 2D: change in direction over time 
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Figure 19: Direction of movement over time for the same simulation. 
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Figure 20: Plots of escape position for 200 uniformly spaced choices of circle radius from 0.1 - 20, 100 
independent simulations per radius. 
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Figure 21: Histograms comparing the distribution of exit angles for radii r = 20, 50, 100 and 200, 
based on 100,000 simulations each. 
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Figure 22: Mean escape times plotted against radius of the circle, based on 10,000 simulations. 
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Quantity 


Description 


Value 




Decay rate of activator 


2 X 10-^ 


n 


Decay rate of global inhibitor 


3 X 10-2 


Tc 


Decay rate of local inhibitor 


1.3 X 10-2 


Da 


Diffusion Coefficient of activator 


4 X lO-'^ 


Db 


Diffusion Coefficient of global inhibitor 


4 


Dc 


Diffusion Coefficient of local inhibitor 


2.8 X 10-6 


ha 


Basal activator production rate 


1 X 10"^ 


he 


Local inhibitor production rate 


5 X 10-3 


Sa 


Saturation CoefRcient 


5 X lO-'' 


•S,. 


Michaclis-Ahnil cu const aul 


2 X 10""^ 



Table 1: Table of Reaction-Diffusion equation constants 
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